function [x] = bandthree(n,a,b,c,y)
% band matrix inversion for 3 band matrix
% copied to matlab from Wiki's fortran routine

%Inputs:
% n - dimension of matrix
% a,b,c - band matrix coefficients
% y is the known variable and x is the variable to solve
% z is dummy var

% matrix deconvolution
Bi(1) = b(1);
Ci(1) = c(1)/Bi(1);
for i = 2:1:n
    Ai(i) = a(i);
    Bi(i) = b(i) - Ai(i)*Ci(i-1);
    Ci(i) = c(i)/Bi(i);
end

% forward substitution
z(1) = y(1)/Bi(1);
for i = 2:1:n
    z(i) = (y(i)-Ai(i)*z(i-1))/Bi(i);
end

% back substitution
x(n) = z(n);
for i = n-1:-1:1
    x(i) = z(i)-Ci(i)*x(i+1);
end

end